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ABSTRACT 

Context. Massive stars influence their environment via stellar winds, ionising radiation and supernova explosions. This is signified 
by observed interstellar bubbles. Such "feedback" is an important factor for galaxy evolution theory and galactic wind models. The 
efficiency of the energy injection into the interstellar medium via bubbles and superbubbles is uncertain, and is usually treated as a 
free parameter for galaxy scale effects. In particular, since many stars are born in groups it is interesting to study the dependence of 
the effective energy injection on the concentration of the stars. 

Aims. We aim to reproduce observations of superbubbles, their relation to the energy injection of the parent stars and to understand 
their effective energy input into the interstellar medium (ISM), as a function of the spatial configuration of the group of parent stars. 
Methods. We study the evolution of isolated and merging interstellar bubbles of three stars (25, 32 and 60 M ) in a homogeneous 
background medium with a density of 10 m p crrT 3 via 3D-hydrodynamic simulations with standard ISM thermodynamics (optically 
thin radiative cooling and photo-electric heating) and time dependent energy and mass input according to stellar evolutionary tracks. 
We vary the position of the three stars relative to each other to compare the energy response for cases of isolated, merging and initially 
cospatial bubbles. 

Results. Due to mainly the Vishniac instability, our simulated bubbles develop thick shells and filamentary internal structures in 
column density. The shell widths reach tens of per cent of the outer bubble radius, which compares favourably to observations. More 
energy is retained in the ISM for more closely packed groups, by up to a factor of three and typically a factor of two for intermediate 
times after the first supernova. Once the superbubble is established, different positions of the contained stars make only a minor 
difference to the energy tracks. For our case of three massive stars, the energy deposition varies only very little for distances up to 
about 30 pc between the stars. Energy injected by supernovae is entirely dissipated in a superbubble on a timescale of about 1 Myr, 
which increases slightly with the superbubble size at the time of the explosion. 

Conclusions. The Vishniac instability may be responsible for the broadening of the shells of interstellar bubbles. Massive star winds 
are significant energetically due to their - in the long run - more efficient, steady energy injection and because they evacuate the space 
around the massive stars. For larger scale simulations, the feedback effect of close groups of stars or clusters may be subsumed into 
one effective energy input with insignificant loss of energy accuracy. 

Key words. Galaxies: ISM - ISM: bubbles - ISM: structure - hydrodynamics - Instabilities 



^ | 1. Introduction lems for galaxy ev olution models (lElmegreen & BurkertH20l61 

. , , • ox, \ . [Scannaoieco et al. 20TI]): details of the feedback prescription 

Feedback from massive stars and their supernovae (SNe) is ch ^ ^ masf , fe a factQr of a few and ±e starg 

^ ■ commonly agreed to be an important driver of the evolution m tQQ concentrated for all prescriptions, which leads to un- 

of galaxies, notab y disk galaxies such as the Milk y Way (e.g. realisticall dec lining rotation curves. In fact, even the mor- 

I Scannapieco et alj2008t | Piontek & Steinmetz] | 201 1|) . It remains hol ical t of a al i e whether or not it has a stel . 

unclear if and under which circumstances such feedback is pos- ^ seems tQ d d Qn ^ feedback implementation. 

itive, i.e. enhances the star forming activ ity, or negat ive^. ft ^ ^ ^ such feedback ^ . ^ for ^ ^ 

decreases or halts star forming activity (|Powell et al. | |2011l. , me 

most-massive acti ve galaxies, where supermassive 

\ylrtra/M7af* 1 1 vamn mnn unnlanr I-i/-\itt nn/»n TaaHrtnnl^ rtnnnrr <" \ «"\ i t-\ ^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 



Moreover, it remains unclear how such feedback occurs on in- black holes m ay domin ate (e. g .lKrausel2005HCroton et alJ2 006 

termediate (temporal and spatial) scales, i.e. in between the i NeS vadba et al l l2008t iKrause & Gaiblerf 1201(1 iGasr^riettl 

immediate surroundings of individual stars - with their radi- 201 9 ' Silk & Mamorf20f2l) 
ation and wind impacts or supernovae - and the large galac- 
tic scales of kpc. Also it is an open question if giant molec- Unless the stellar ejecta are permanently removed from their 

ular cloud lifetimes of the order of 0.1 Gyr can be recon- host gal axy, which may be pa rticularly relevant for smaller 

ciled with this feedback. These uncertainties cause severe prob- galaxies (|Lanfranchi et_aL||2006|), chemical evolution in general 

is not very sensitive to these details of the feedback descrip- 

* E-mail: Martin.Krause@universe-cluster.de tion: The models reproduce observations well with an instan- 
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Fig. 1. Comparison of CU.CT (solid line) and the HLLD.CT 
(dashed) methods for the blastwave test problem. The problem 
is solved on a 3D Cartesian mesh in each case, using two AMR 
levels. The common snapshot time is 0.0012. The contact sur- 
face (at «20/60 on the horizontal axis) is essentially identical. 
For the backwards shock (at s;36/44) CLLCT leads (i.e. is closer 
to the centre) by about one cell. CLLCT also leads for the for- 
ward shock (at « 12/68), but only by a fraction of a cell. 



taneous mixing approximation (e.g. iPraritzo s & Silk 1998]) and 
even if one considers mixing delays due to having the stellar 
ejecta in a hot pha se or in galactic fountains for some time 
dSpitoni et al.ll2.Q09b . Yet, such details do matter in the context 
of prompt metal enrichment by ne arby massive stars, as dis- 
cussed for the protosolar nebula (e.; 



. .Gounelle & Meiboml2 008; 
2011b . Also, in order to un- 



iGounelle et afll2009l: iMakide et all 
derstand the spread in the observed ratio of r-process (produced 
in core collapse supernovae, only) to ^-process elements (pro- 
duced later by less massive stars) dSimmerer et al.ll2004l) . and 
especially the /--rich metal poor stars dSneden et al.ll2008l) . feed- 
back details might be important. 

The properties of the interstellar medium (ISM), i.e. 
its morphology wit h imprinted b u bbles and superbubbles 
dGruendl et al.l l200Ct lArfhurl 120071: Ichul 120081: ISasaki et al.l 
1201 ll) as well as molecular-cloud fragments in formation or 
in dispersal, and its level of turbulence, are strongly af- 
fected by the physics and dynamics of stellar feedback (e. g. 



de Avillez & Breitschwerdtl 12004 120051 iDobbs etaP 1201 lblat 
Ntormousi et alJuOl lb . The actual agents of stellar feedback 
are massive stars, born in the dens er parts of the interstel- 
lar medium (for recent r eviews see iMcKee & Ostrikerl 120071: 
IZinnecker & Yorkell2007l) . The interaction via winds and ion- 
ising radiation of a single massive star wit h its surroundings i s 
usually referred to as "interstellar bubble" dWeaver et al.lll977l) : 
Strong winds are shocked close to the star and produce a hot 
overpressured bubble, which drives an expanding shell of swept- 
up, shocked ambient gas. The shell may be partially or com- 
pletely ionised by the Ultraviolet emission of the central star. 
The increased pressure may additionally push the leading shock 
front. Interstellar bubbles are usually not energy conserving, be- 
caus e of the radiative l osses of the shocked ambient medium 
(e.g. lWeaver et al.lll977 | ). This has been nicely demonstrated by 
the observations o flGruendl et al.l (12000): For some bubbles, they 
resolve the radiative leading shock wave, with the highly excited 
[ O III ] tracing the hottest outermost gas, and Ha tracing a some- 




60 Ms 



32 Ms 



25 Ms 



10 



Time / Myr 



Cumulative mass input 




a 30 r 



4 6 
Time / Myr 



10 



Fig. 2. Energy (top) and mass (bottom) input. We use the output 
of a 25 M (dashed, labels: Ms= M ), 32 M Q (dotted) and a 
60 M Q star as input for our simulations, separately or combined. 
The difference between the total mass output and the initial mass 
indicates the mass of the dark remnant. The energy is given in 
"Bethe" = 10 51 erg. 



what cooler surface inside of [O III]. This indicates that the 
leading shock front in these cases is shock ionised rather than 
photo-ionised. Radiative energy losses are substantial, but hard 
to quantify in detail (e.g. iGarcia-Segura & Mac Lowlll995t) and 
affect wind-blown and supernova related bubbles alike. 

Many molecular clouds host massive stars in groups. The 
bubbles of these stars have to interact, because the sizes 
of the individual b ubbles (parsecs, e.g. IWeaver et al.l 119771 : 
iGruendl et al.l l2000h i s comparable to the si ze of the par- 
ent molecular clouds dKainulainen et alJ 1201 ll) . Also, in star 
forming regions, smaller groups of stars are oft en located 
withi n distances of tens of parsecs (e.g. Orion, IVoss et all 
120101) . The interaction of individual bubbles leads to the for- 
mation of superbubble s (|Tenorio-Tagle & Bodenheimerl 119881 : 
lOev. Clarke. & Massevl 120011: IChul 120081: lOevI l2009l for re- 
views): The expansion of the combined superbubble is of- 
ten described by the same model as for individual bubbles, 
which predicts expansion rate and shell size, if the energy in- 
put and the ambient density are k nown. Superbubbles may reach 
sizes of hundreds of parsecs (e.g. | Tenorio-Tagle & Bodenheimerl 
ll988HBreitschwerdt & de Av illez 2006t: ISasaki et al.ll201 ll) . But 
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Table 1. Simulation parameters 



Label 


Star mass /M 


X/pc 


Y/pc 


Z/pc 


Res. / pc 


/Jo/cm 3 


T /K 


S25 


25 











2.1 


10 


121 


S32 


32 











2.1 


10 


121 


S60 


60 











2.1 


10 


121 


3S0 


25 











2.1 


10 


121 




32 



















60 

















3S1 


25 


-30 


10 


10 


2.1 


10 


121 




32 
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-10 













60 

















3S2 


25 


-60 


20 


10 


2.1 


10 


121 




32 


50 


-10 













60 

















3Sl-mr 


25 


-30 


10 


10 


1.0 


10 


121 




32 


-25 


10 













60 

















3Sl-hr 


25 


-30 


10 


10 


0.52 


10 


121 




32 


-25 


10 













60 


















they often appear to be too s mall and too bright in X- 
rays compared to models (e.g. lOev & Garcia- Segural EZOOH 
iJaskot et all 1201 lb . Possible explanations include energy dis- 
sipation due to mass loading or uncertainties in the stellar 
wind data, e.g. due to clumping. There is considerable dis- 
agreement also in the litera ture about the deta ils of the physi- 
cal model of superbubbles. IJaskot et"a"D d201 lh argue for mass 
evaporation by thermal conduction for the X-ray bright super- 
bubbles DEM L50 (N186) and DEM LI 52 (N44) in the Large 
Magellanic Cloud (LMC), because the X-ray temperatures are 
low and there are no direct detections of clouds being ablated. 
Yet. IBreitschwerdt & de Avillezl (1200 61) present a detailed hydro- 
dynamic simulation without thermal conduction, which repro- 
duces much of the available observational data of the local bub- 
ble, which comprises the Solar System, and the nearby Loop 1 
superbub ble. Most notably the O VI absor ption column pre- 
dicted by IBreitschwerdt & de AvilleJ (l2006h, whic h is strongly 
produced by evaporation (e.g. Weaver et al.lll977l) agrees well 
with observations. 

Understanding of the physics of bubbles and superbubbles 
is the key ingredient in order to gauge the efficiency of stel- 
lar feedback. It is of particular importance to assess the effec- 
tive energy input into the ISM. Our group has embarked on this 
task, and has synthesised the total energy input into molecular 
clouds for realisti c stellar populat ions based on recent stellar 
evolution models dVoss et al.ll2009l) : Averaged over all massive 
stars (8M Q < M < 120M o ), the energy input due to winds is of 
order 10 50 erg/star. Supernovae contribute about ten times more. 
The energy injection is extended over several tens of Myr and 
has a peak near four Myr with a shallow decline afterwards. 
Winds dominate before the peak and supernovae afterwards. 
Substantial variations from cluster to cluster are expected due 
to the sparser sampling at the massive end of the initial mass 
function. 

Stars are born in the densest regions of the ISM. Much of 
the injected energy is therefore quickly lost to radiation in cool- 
ing shock compressed shells. Hydrodynamic simulations have 
been used to assess the effective energy input into the ISM. The 
energy deposition efficiency of isolated massive stars in their 
wi nd phase has been assesse d in 2D hydrodynamic simulations 
bv iFrever etall d2003l [2006b . Though they include the effect of 
photo-ionisation, they show that the gas dynamical effects are 



dominated by the mechanical energy input: For example, for a 
35 (60) M star, they expect only 17 (5) per cent of the energy 
transfered to the ISM in their simulation being due to the ef- 
fect of ionisation. They give their energy deposition efficiency 
as fractions of the radiative energy input. Scaled to the mechani- 
cal energy input, they find that about 38 (9) per cent of the input 
energy has been added to their ISM at the end of their simula- 
tions for the 35 (6O)M star. The dynamics of two wind bubbles 
(25 M a and 40 M stars) sep arated by 16.2 pc has been stud- 
ied bv lvan Marie et all (120121) . The two bubbles quickly merge, 
sweeping the colliding parts of the wind shells away into the 
bubble of the lower mass star, due to the pressure difference in 
the bubbles. An aspherical superbubble is then formed, which 
isotropizes after a few Myrs. More interesting details are ob- 
served which we refer to below, when we c ompar e them with 
our findings in Section [3] INtormousi et al.1 (1201 ll) have simu- 
lated the merging of two superbubbles in 2D with identical stel- 
lar cont ent. One of the most int eresti ng findings in the simu la- 
tions of INtormousi et all (1201 ll) and Ivan Marie et al.l (1201 2l) is 
the o ccurrence of the Vishniac thin shell instability (IVishniaci 
119831). This ins t ability is stro ngly suppressed in the simulations 



ty is stro 

of lFrever etlill (l2003l 120061) due to the thickening of the shell 
because of the increased pressure due to the ionisation. The 
Vishniac instability is interesting as it may create observable fila- 
mentary features, and thick filamentary shells, and thus discrim- 
inate between models. 

Here, we address the effective energy injection into a homo- 
geneous ISM for three interacting interstellar bubbles with 3D 
hydrodynamics simulations, using standard ISM thermodynam- 
ics. We neglect the effect of ionising radiat ion, because i t is ex - 
pected t o be less important in this context dFrever et al.l d2003l) . 
see e.g. iGritschneder et all d2010l) for the effects of ionising ra- 
diation). We take as our starting element a group of three coeval 
massive stars, 25, 32 and 60 M Q , respectively. We study the de- 
position efficiency of energy injection as a function of distance 
between the stars. We find a hig h efficiency in the wind p hase, 
comparable to the 2D re sults of Frever et al. j|2003l 120061) and 
the 2D and 3D results of Fierlinger et al.l d2012l) . details of bub- 
ble merging similar to Ivan Marie et al. "(12012 ) and an enhance- 
ment of the feedback efficiency by about a factor of two for 
grouping of the stars closer than about a few tens of pc. The 
energy of supernovae that explode within superbubbles is dis- 
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Fig. 3. Column density integrated over Z-direction and ^-direction (left and middle columns, respectively) and midplane density 
(right column) for three different snapshot times from top to bottom for run 3Sl-hr. The projections of the three massive stars into 
the X-Y plane is indicated as small red stars in the density plots on the right. The 60 M Q star blows the biggest bubble from the 
origin. The 32 M bubble towards its lower left (XK-plots) is only slightly bigger than the one of the 25 M star above. The shell 
forms spikes and dense clumps due to the combined action of Vishniac and thermal instability. A movie is provided with the online 
version. 



sipated on a timescale of about 1 Myr. Additionally, we show 
the first column density renderings of prominently Vishniac un- 
stable 3D shells, which should give a first approximation of the 
observational appearance of the Vishniac instability. 



2. Simulations 



We carry out 3D hydrodynamic simulations with the 
Nirvana 3.5 code (IZiegrer||2008l l201lb . evolving the conserva- 
tion equations for mass, momentum and energy. Nirvana 3.5 is a 
conservative, finite volume code and combines block structured 
adaptive mesh refinement (AMR) with parallelisation by the 
message passing interface (MPI) library. 
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Fig. 4. Figure[3]continued, but with all scales adapted to the snapshots presented in this figure. 



2. 1 . Numerics and code tests 

The main solver modules are an HLLD solver (HLLD_CT), 
applying the I D app roximate Riemann solver of 
iMivoshi & Kusanol (I2005I) dimension-by-dimension in 3D, 
and a se cond-order Ce ntral-Upwind scheme (CLLCT, full 
details in IZieglerl |201 ll) . We work in Cartesian coordinates 
throughout. In order to check the isotropy of the solution in this 
geometry and also for differences between these solvers, we 
have re-run and analysed the adiabatic blastwave test problem 
that comes with the code with both solvers (Figure [TJ. Here, a 
fixed amount of thermal energy is initially deposited in a finite 
circular region of 22 cells diameter. In both cases, a reasonably 
spherically symmetric bubble develops, with a forward shock, 



a contact surface and a backward shock. The contact surface 
evolves identically for both solvers, forward and backward 
shock are lead by the solution of the CLLCT solver by at most 
one cell. Hence, both methods yield a very similar result for 
symmetrical bubble expansion. For the same solver but different 
angular directions other than the grid axes, the radii of the 
different features of interest differ by typically one and up to 
about three grid cells. 

We have initially selected the HLLD XT solver but encoun- 
tered severe vacuum formation problems (very low pressure) 
near contact surfaces for our high resolution runs. For all the 
simulations presented in this article, we have therefore employed 
CLLCT. 
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We use standard ISM thermod ynamics with r adiativ e cool- 
ing and photo-electric heating (see iPiontek et al.l (120091) for de- 
tails), employing the standard iterative procedure of Nirvana 3.5. 
Cooling is always strong for our wind shells, which tend to get 
thin and eventually also Rayleigh-Taylor and Vishniac-unstable. 
The instabilities evolve differently for different flux limiters: 
Test simulations with all flux limiters provided (minmod, super- 
bee, monotonised-centred, and Van Leer) showed that for the 
monotonised-centred and the superbee limiters, the instabilities 
are systematically different for parts of the shell which move par- 
allel and diagonal to the grid axis. Van Leer and minmod both 
yield almost isotropic results at our highest resolution, at the ex- 
pense of being more diffusive, as expected. We have correspond- 
ingly adopted the minmod flux limiter. Nirvana 3.5 offers an ad- 
ditional multi-dimensional limiter which we also use, and where 
we have adjusted the parameter experimentally to yield optimal 
isotropy for shell instabilities. 



2.2. Setup 

The computational domain is a cubic Cartesian grid, 400 pc on a 
side resolved by 24 cells for the base level. The mesh is refined 
whenever a combined threshold of first and second derivative for 
density or respectively velocity is exceeded. Additionally, we al- 
ways keep the wind injection region at the highest refinement 
level. Effectively, the wind shell and everything inside is always 
refined to the highest level. For most of our runs we use three 
levels of adaptive mesh refinement, which would correspond to a 
uniform grid of 192 3 cells with a resolution of 2. 1 pc. Simulation 
3S 1 -mr and 3S 1 -hr use four and five levels of refinement, result- 
ing in 1 and 0.5 pc resolution, respectively. Boundary conditions 
are formally periodic, but we only use data from snapshots where 
the shells are entirely contained in the computational domain. 

We fill the grid initially with a homogeneous medium. Then 
we choose one (three) injection regions of eight pc radius in ev- 
ery case. Each injection region gets assigned a star of a partic- 
ular mass. We inject mass and thermal ener gy according to the 
stellar evolutionary tracks of rotating stars of | Mevnet & Maederl 
( 2005b and wind veloci t ies fr om lLamers et al.l ( 119951) and 



Evolution of maximum density 



iNiedzielski & Skorzvnski d2002l) for the Wolf-Rayet phase, as 
compiled in lVoss et al.1 d2009t) . We use 25, 32, and 6OM stars, 
with supernovae at 8.6, 7.0 and 4.6 Myr, respectively. The time 
resolution of the stellar evolution table is 0. 1 Myr. Mass and en- 
ergy input are shown in Figure [2] The mass density is initially 
set to 10 m p cm -3 everywhere in the computational domain. The 
temperature is set in equilibrium between cooling and heating, 
121 K. All velocities are initially zero. More details for each in- 
dividual run are provided in Table Q] 

3. Results 

The time evolution of our high resolution run 3Sl-hr with three 
stars at different locations is shown in Figures|3]and|4] At a given 
time, the bubble size increases monotonically with the mass of 
the parent star, with the central 60 M G bubble dominating the 
gas dynamics. As expected, the shocked ambient medium cools 
very quickly and consequently gets compressed into a thin shell 
for each bubble. The shell is subject to a combination of thermal 
and IVishniad (1 19831) instabilities^. The bubbles start to merge at 
around 2 Myr. At the first snapshot in Figure |3]( 1-95 Myr), the 



1 Although we have carefully chosen the flux limiter, the shell insta- 
bility evo lves still somewhat ani sotropically. This is similar to the 2D 
results of iNtormousi et all d201 if) with the RAMSES code, where even 
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Fig. 5. Maximum density as a function of time for runs 3S1- 
hr (solid), 3Sl-mr (dashed) and 3S1 (dotted). The horizontal 
lines correspond to the critical compression above which the 
Vishniac instability is triggered for supernova (lower line) and 
wind (thicker upper line) shells according to IVishniac & Rvul 
(I1989I) . The axis on the right shows the overdensity factor over 
the undisturbed ambient medium. 

3S1 -hr: slice at Y,Z = 1 pc, T = 8.53 Myr 




Fig. 6. One-dimensional slices in X-direction through run 3S 1 -hr 
at time T = 8.53 Myr. The Y and Z coordinates are chosen appro- 
priately for the slices to include the position of the only remain- 
ing star at that time (25 M Q , at X = -30 pc, indicated by the star 
in the middle diagram). Top: positive X-velocity (blue, dashed 
line), negative x-velocity (red dash-dotted line) and sound speed 
(solid black). Middle: pressure (logarithmic) . Bottom: density 
(logarithmic). See text for details. 



shell interface between the 60 M bubble and the 32 M Q bubble 
has just burst. Up to this point, each bubble has had its individ- 
ual bubble pressure, which is largest for the 60 M Q bubble. Its 
hot gas can be seen to stream through the hole in the shell. The 
shell in terface then behave s much like a cloud, being ablated by 
a wind dPittard et al.ll2005l) : Kelvin-Helmholtz instabilities at the 
contact surface lead to mixing of the cloud gas into the hot phase. 
The shell interface has completely dispersed until the next snap- 
shot at 4.05 Myr. We have checked the effect of different flux 



a five times higher spatial resolution could not get the shell instabilities 
completely isotropic. 
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limiters in this phase: Less diffusive ones allow smaller holes, 
which delays the erosion process compared to the more diffu- 
sive case. The final results are however very similar. 

The density slice at 4.05 Myr shows the weaker winds of 
the smaller stars to be pushed aside by the one of the most mas- 
sive star. The larger part of the 60 M Q bubble remains unaffected 
by the action of the smaller stars. The 60 M Q star explodes at 
4.6 Myr. The sudden energy injection due to the supernova com- 
presses the shell further (Figure [5]) and accelerates it, triggering 
the Rayleigh-Taylor instability (RTI). The RTI may cause fila- 
mentary structure inside the shell. Also, the outwards directed 
flow field, centred around the most massive star before its ex- 
plosion, is no longer present. Thus, from this time on, we find 
filamentary gas inside the shell, seen in the individual density 
slices. The effect of the winds of the smaller stars in this phase 
can hardly be noticed. The second supernova (7.0 Myr) leads 
to a further acceleration and compression of the shell, causing 
more RTI filaments. The snapshot at 8.53 Myr shows the super- 
bubble when 2 stars have exploded already, and the third is in its 
Wolf-Rayet phase. This snapshot demonstrates nicely that our 
ansatz with thermal energy injection may also cope with situa- 
tions when the backward shock within the stellar ejecta is un- 
usually far from the star: One can clearly see the declining den- 
sity away from the star due to adiabatic expansion (ID-slices in 
Figure [6]). The wind turns supersonic immediately outside the 
driver and shocks roughly 20 pc away from the star. A second 
structure is visible at varying distance from the star, up to about 
50 pc: This is what we would expect to be the forward shock in 
the standard picture. Due to the high ambient pressure, it is only 
a sound wave. The pressure inside of this structure is slightly 
reduced due to the ongoing expansion. The final supernova at 
8.6 Myr causes again mass entrainment into the bubble due to 
the RTI. The bubble then keeps expanding with decreasing inte- 
rior density fluctuations until the end of the simulation at 15 Myr. 

The highest densities in the shell, around 180 times the ambi- 
ent density, are reached for roughly 1 Myr after each supernova, 
where for the later two supernovae, the compression peaks have 
merged (Figure|5]). At late times the density increases again (see 
below for details). We show a zoom on the highest density re- 
gion in the final snapshot in Figure [7j The density maximum is 
located in the dense shell, where two humps of the Vishniac in- 
stability (compare Section [37X1 below) cross, and more towards 
the interior of the bubble. The velocity field in the shell is still 
dominantly outwards with substantial Mach numbers. Yet, prob- 
ably enhanced by the large scale vortices which dominate the 
shell interior at that time, there is also some non-radial mo- 
tion. The slightly converging velocity field has to be responsible 
for the high density, as the region is substantially overpressured 
compared to the environment. At earlier times (compare above), 
such maxima in density and pressure could have been in pressure 
equilibrium with their surroundings. At this late time, the bubble 
interior is already underpressured with respect to the environ- 
ment, and so we expect that the maximum is temporary, unless 
such clumps become self-gravitating. This seems quite likely, 
given the pc-scale size, low temperature (below 20 K) and high 
mass (few hundred M Q ) of the clump (Jeans length: « 2 pc). 
Yet, self-gravity is not included in the simulations and therefore 
details, such as triggered star formation, are beyond the scope of 
this article. 

3.1. Vishniac instability 

The shells are subject to various instabilities. The Rayleigh- 
Taylor instability is especially prominent during the strong ac- 
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Fig. 7. Shell details for the final snapshot of run 3Sl-hr. Shown 
is an X-Y zoom of density, pressure, temperature and Mach 
number, as indicated on the individual panels, around the po- 
sition of the maximum density, which is located at (X, Y,Z) = 
(-30,-79,-52) pc. Velocity vectors are overlaid on the density 
plot. The high density region is overpressured and has a temper- 
ature below 20 K. See text for more details. 



celeration phases after each supernov a. Because this is the first 
time that the Vishniac shell instability dVishniadll983l) is promi- 
nently seen in 3D simulations, we give a few more details. The 
Vishniac instability is an overstability: differences in column 
density for adjacent regions of a shell cause gas flow from the 
high column density region into the region with smaller column 
density. This continues in general until the situation is reversed 
and the region with initially smaller c olumn density finally has 
the greater one. IVishniac&Rviilll 1989ft derive a critical overden- 
sity for the shell over the unshocked ambient gas of a factor 10 
and 25 for a blastwave with initial energy injection and constant 
energy injection rate, respectively, to become unstable, such that 
the peak density increases in eac h cycle. The shell then devel- 
ops a characteristic spiky pattern dNtormousi et al.|[2.Ql U iDrakd 
1201 2l Figure[3}, in density slices. The 3D structure of the shell is 
granular with a regular filamentary pattern (Figure[3]and[4]). The 
regularity is of course related to the grid structure, because this 
is the most important perturbation. In column density, we find a 
web formed of polygons. These polygons have typically four to 
six sides. The sides are however not always aligned with the co- 
ordinate axis or the diagonals, and some are clearly curved. The 
typical polygon diameter is about 10 pc. At the intersections of 
the filaments, density and column density achieve their highest 
values. These points lag behind the shell. Particularly high den- 
sities may be achieved, when left and right part of an inwards 
spike merge. This seems to have happened for the density max- 
imum at the final snapshot we show in Figure |7j But from a 
detailed inspection of several snapshots, we conclude that this 
should happen frequently. 

We show the peak density over time in Figure [5] Clearly, 
the densest parts of the shell of run 3Sl-hr satisfy the criteria of 
IVishniac&Rviilll 1989ft from before 2 Myr throughout the simu- 
lation, in agreement with Figure[3] The low reso lution simulation 
3S1 ge nerally stays below the wind criterion of lVishniac & Rvul 
(119891) . Corresponding ly, the Vishniac instability is much less 
pronounced (Figure[8). Mac Low & Normanl (1 19931) have shown 



7 



Martin Krause et al.: Emergence of superbubbles - energy efficiency & Vishniac instabilities 



100 



CL 



3S1, time: 12.69 Myr 



Hi 



3S1-hr: mass-weighted density vs. nan-radial Mach number histograms 



-100 



■1 ,34 



■1 .51 



■1 .69 



1.86 



-100 -50 50 100 
X/pc 

3S1-hr, time: 12.56 Myr 



100 



CL 




-100 



1 .25 



■1 .48 



■1 .72 



1.96 



E 
a 
m 



in 
c 
v 
■a 
c 
G 

3 



□ 



E 

o. 

CT-' 



en 
t 

"D 
£Z 

E 

O 

en 
o 



-100 -50 50 100 
X/pc 

Fig. 8. Column density at a comparable late evolution time for 
runs 3S1 (top) and 3Sl-hr (bottom). The high resolution bubble 
is more spherical, larger, achieves higher peak column densities 
and the Vishniac instability is more pronounced. 



that the instability is connected to transonic motions in the shell 
perpendicular to the expansion direction. We evaluate these non- 
radial velocities for the undisturbed (with respect to the interac- 
tion of the bubbles of the other two stars, here we use X > 0) part 
of the shell in Figure [9] The 2D mass weighted histogram over 
logarithmic density and non-radial Mach number, with respect 
to the local speed of sound, shows that only dense shell gas ac- 
quires substantial non-radial Mach numbers. At high densities, 
indeed most of the gas has Mach numbers around and below 
unity. 
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Fig. 9. Analysis of the non-radial Mach number, i.e. the Mach 
number perpendicular to the direction of the shell's expansion. 
Only the region with positive X-coordinate, which corresponds 
to the undisturbed part of the 60 M© bubble, is taken into ac- 
count. Left: 4.05 Myr, right: 7.47 Myr. The upper parts show 
the non-radial Mach number versus the logarithm of the den- 
sity. Colour encodes the mass per bin, where each bin spans 
0.05 in Mach number and 0.06 dex in logarithmic density. No 
appreciable non-radial motions are found for the hot bubble in- 
terior, whereas the dense shell material shows Mach numbers 
of order unity. The lower parts show mass weighted non-radial 
Mach number histograms (vertically collapsed versions of the 
plots above). The plots are dominated by the quiescent ambient 
medium. The mass with given non-radial Mach number declines 
strongly around a Mach number of unity towards higher Mach 
numbers as expected for shells dominated by the Vishniac insta- 
bility. 



tion. We define the response H to be the energy retained in the 
ISM divided by the input energy: 



£ism(0 - £ism,o 



(1) 



where we define as ISM the whole gas present in the computa- 
tional domain, including the hot bubble interiors with their stel- 
lar ejecta. 

H is generally of order ten per cent. It is higher whenever the 
energy input rate increases. This is especially well visible at the 
time of the three supernovae at 4.6, 7.0 and 8.6 Myr. Here, fl 
reaches peak values between 20 to 40 per cent. % is smaller for 
phases of decreasing energy input rate. This is particularly well 
visible after a supernova. About 1 Myr after each supernova, % 
drops to roughly five per cent. The characteristic decay time of 
the retained energy increases for each consecutive supernova. 
When the energy input ceases, the ISM energy is lost to radiation 
on a timescale of Myrs, with fi dropping to 2 per cent roughly 
4 Myrs after the last supernova. 

Steady, continuous energy injection is clearly more effective 
in energising the ISM than sudden bursts such as from infrequent 
supernovae. 



3.2. Energy evolution: general observations 

We show the total input energy over time together with the en- 
ergy retained in the ISM where the initial thermal energy is sub- 
tracted in Figure [10] The retained energy is generally below the 
input energy because the gas is initially in radiative equilibrium 
and suffers net radiative losses during the course of the simula- 



3.3. Resolution effects 

We have repeated run 3Sl-hr at a half and a quarter of the 
original spatial resolution. Morphologically, the bubbles are less 
spherical, smaller and the Vishniac instability is less developed 
at lower resolution (Figure |H). We compare the energy evolu- 
tion of the three runs in Figure QT| The retained energy differs 
by much less than a factor of two between simulations at dif- 
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3S1-hr: input and retained energy 
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Fig. 10. Top part: input (solid) and retained (dotted) energy for 
run 3Sl-hr. The response 7? (retained energy divided by input 
energy) is shown in the bottom part. See text for more details. 



3S1, 3S1-mr&3S1-hr: energy tracks 
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Fig. 11. Resolution effects on the retained energy. Top part: re- 
tained energy for run 3Sl-hr (solid line, high resolution), 3S-mr 
(dashed line, intermediate resolution) and run 3S1 (dotted line, 
low resolution). The bottom part shows the energy ratio 3S1- 
hr/3S 1 -mr (solid line), 3S 1 -mr/3S 1 (dashed line) and 3S l-hr/3S 1 
(dash-dotted line). In each case, the data for the higher resolution 
run has been interpolated to the data output times of the lower 
resolution run. The spikes at the supernova times are artefacts of 
the interpolation process at the discontinuities of the functions. 
The horizontal dashed line indicates equality for comparison. 
The energy increases similarly for each doubling of resolution. 
The general functional behaviour is independent of resolution. 
See text for more details. 



diated away. Also, the peak density at a given time depends 
strongly on resolution (Figure [5}, which also changes the ther- 
modynamics. 

3.4. Energy evolution: varying stellar distances 

We have carried out a set of simulations, where we varied the 
positions and distances of the same three stars (Figure H2t . 
Because of computational limitations, these simulations have 
been carried out at 2.1 pc resolution. This is physically justi- 
fied by the convergence of the general shape of the energy tracks 
(Figure [TT). For obtaining the large distance limiting case, we 
have simulated each of the three bubbles in a separate simula- 
tion (S25, S32 and S60), and added their energy tracks for com- 
parison to the other cases. We model the closely-spaced extreme 
case, where the bubbles have merged instantaneously, by putting 
the driver regions of the three stars on top of each other at the 
grid origin (3S0). Additionally, we performed two simulations 
with intermediate star positions (compare Table [TJ, where we 
actually observe the bubble merging during the simulations (3S1 
and 3S2). 

We see small differences in the energy tracks during the first 
=s 0.5 Myr. They are expected because during this time, the driver 
region is evacuated and the bubble shape is established. It makes 
of course a difference, if the three stars share the same driver re- 
gion (3S0), or if each star has its own. Also shifting the driver 
region on the grid makes the volume of the individual driver re- 
gions slightly different, by a few per cent, due to resolution ef- 
fects at the driver boundary. This translates to a few per cent 
difference in total energy, which is visible in Figure[T2l(bottom). 

Once the bubbles are established properly on the grid, i.e. 
after about 0.5 Myr, all configurations have essentially the same 
energy response until the first supernova at 4.6 Myr. The reason 
for this is the predominance of the energy injection of the 60 M & 
star. The energy tracks begin to differ slightly after the first star 
has exploded. The divergence increases abruptly after each su- 
pernova. But for very long times after the final explosion, the 
tracks converge again towards a common value. 

Among the four configurations, the energy varies at times by 
up to a factor of three. A typical value after the second supernova 
is a factor of two. Throughout the simulation time, the energy 
is essentially highest for run 3S0 (all stars at same place) and 
lowest for very large distance (sum of S25, S32 and S60). The 
two configurations with intermediate distances, where the bub- 
bles merge during the respective simulations, show intermediate 
energies. The run where the bubbles merge early (3S1) behaves 
almost identical to the case where the driver regions are on top 
of each other (3S0). 



ferent resolution. The differences are more pronounced at later 
simulation times. Finer spatial resolution always leads to more 
energy in the ISM. For an increase of the resolution by a factor 
of two, we find an increase of the retained energy by 20-30 per 
cent. This agrees with the greater bubble diameter at higher res- 
olution (Figure|8]l. The overall functional behaviour is very well 
converged. 

The reason for the changes with resolution is very likely 
the details of the shell evolution. In the absence of other per- 
turbations, instabilities are triggered on the resolution level. 
Additionally, the Vishniac instability is only marginally devel- 
oped at low resolution. This might lead to more non-radial ki- 
netic energy at higher resolution, which is not immediately ra- 



3.5. Shell widths 

We find that our simulated shells are widened due to the Vishniac 
instability. For the determination of the shell width, we average 
the column density maps over the angle, and identify the shell 
as radial interval where the column density is at least five per 
cent higher than in the undisturbed medium. The shell width is 
shown in Figure Qj] as a function of time and radius, respec- 
tively, for runs 3S1 and 3Sl-hr. For this analysis, we only use 
late snapshots, where the superbubbles are well established. 

The shell width is typically in the tens of per cent regime 
and increases with time. The result does not depend on the reso- 
lution. 
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3S1-hr: Averaged column density, time: 12.56 Myr 



Energy tracks: stellar distance variation 




Fig. 12. Energy tracks for different simulations, where only the 
positions of the three stars differ. Run labels in the legends are 
explained in Table Q] S25+S32+S60 refers to the sum of the 
energy tracks of the three simulations of the bubbles of the 
isolated 25 M , 32 M and 60 M stars, respectively, which 
corresponds to a very large distance. 3S0 is the opposite case, 
where three stars are in the same region. Top: Absolute values. 
Middle: Three-stars simulations relative to the sum of the three 
isolated bubbles. Interpolations always use the 3S0 time base. 
Interpolation artefacts are visible at the discontinuities due to 
the supernovae (4.6, 7.0 and 8.6 Myr). Bottom: Difference of the 
very similar energy tracks of runs 3S1 and 3S0, normalised to 
3S0 as a percentage. The solid red line marks zero. The differ- 
ence has been set to zero for the time intervals 10,000 yrs around 
each supernova in order to mask the interpolation artefacts. See 
text for details. 



4. Discussion 

We have investigated the environmental impact of a group of 
three massive stars via 3D hydrodynamic simulation. Herein, 
several assumptions and simplifications were necessarily intro- 
duced: 

We have adopted a uniform background density of 
10 m p cm -3 . On scales of ten pc and smaller, the den- 
sity will in reality be at least a factor of ten higher (e.g. 
iKainulainen et alj201 ll) . On scales of 100 pc, the density should 

(e.g. 



become equal to or even smaller than about 1 m p cm -3 
Ide Avillez & Breitschwe rdt 2005). Hence, our choice should be 
realistic for the tens of pc scales we sim ulate (compare also 
iFrever et alj|2003l 120061 tvan Marie et al.ll2012l) . The real ISM 
has a rich spatial structure, whereas we use a homogeneous dis- 
tribution. This is a significant difference. For a porous ISM, 
the injected wind/SN energy could escape through low density 
regions making the bubbles smaller (Fierlinger et al. 2012b). 
For such a situation, one should also expect pronounced bub- 
ble asymme tries. Indeed, such asym metries are found in obser- 
vations (e.g. lChurchwell et al. 2006|). Yet, in order to be able to 
compare the effect of different spatial configurations of the stars, 
and not to be dominated by local environmental effects it is nec- 
essary to use a homogeneous background density. 
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Fig. 13. Shell width for column density maps. Top: Angle- 
averaged column density over radius for run 3Sl-hr at 
12.56 Myr. From such plots, the shell width has been determined 
as the radial range where the column density is at least five per 
cent greater than at large radii (undisturbed gas). The shell width 
determined in this way is shown in the middle plot as a function 
of time, and in the bottom plot as a function of outer radius. 
Black pluses are for run 3S1, red stars for run 3Sl-hr. The aver- 
age shell width does not depend significantly on resolution. 



Also the thermodynamics has been simplified: we use stan- 
dard ISM heating and cooling functions. In reality, one would 
have to do a full radiative transfer calculation at each timestep, 
as the stars that are responsible for the winds, also have a 
strong photo-ionising output with associated photo-electric heat- 
ing. We use average photo-electric heating of the diffuse UV 
background, but simulate a region close to massive stars, where 
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Fig. 14. Pressure maps of run 3S1 (left) and 3S2 (right), shortly 
after the second (top) and third (bottom) supernova. 



the UV radiation field is much stronger. This should therefore 
result in general in an underestimate of the heating function for 
gas exposed to radiation from a nearby hot star. Similarly, photo- 
ionisation might in reality contribute towards the pressure in the 
shells. Yet, we estimate that the effect on the energetics is small 
(compare SectionQ] above). 

We have found that with the standard ISM thermodynam- 
ics, the peak shell density does not converge with finer reso- 
lution. It is not immediately obvious that this should be so, as 
the photo-electric heating we take into account could in princi- 
ple have produced high enough pressure to limit the shell com- 
pression. Yet, with our highest resolution of 0.5 pc, this has not 
been the case. The shell density results from several effects: At 
a given pressure level, there is a density and a temperature that 
correspond to thermodynamic equilibrium. When the pressure 
in the bubble increases, e.g. because a supernova has happened, 
the shell can however not adjust immediately to the new equilib- 
rium pressure level, because the gas has to be swept together in 
a finite time. The inverse happens at late times (as demonstrated 
in Figure [7]l when the bubble pressure strongly decreases, but 
the clumps in the shell cannot expand fast enough to remain in 
pressure equilibrium. The compression is of course also limited 
by the resolution. The non-convergence therefore means that the 
bubble pressure is high enough for a sufficiently long time so 
that compression of the shell, or some clumps therein, to even 
higher densities may occur if one would repeat the simulation 
with an even higher resolution. For gas on the thermodynamic 
equilibrium curve in the relevant density regime, higher densi- 
ties correspond to lower temperatures. Yet, the temperature in 
the clumps in our simulated shells, frequently reaches tens of 
Kelvin, already, and the coldest clumps reach « 20 K. This cor- 
responds to » 10~ 12 dyn cm" 2 . After a supernova, the bubble 
pressure reaches up to a hundred times more, which would cor- 
respond to sub-K elvin equilibrium tem peratures. Compared to 
observations (e.g. iPreibisch et alJl20l2l) . the ISM in star form- 
ing regions rarely reaches temperatures below about 20 K, and 
20 K to 100 K are typical for the dense phase. Similar temper- 
atures are also found in our simulated shells. Other effects like 
magnetic fields, self-gravity or feedback by the new stars, which 
in reality might form in our dense clouds, may affect the cloud 
compression, but are not included in our simulation. Thus, even 
if the compression would increase still further if one would carry 
out the simulations at yet higher resolution, this would not nec- 
essarily be more realistic, as the high density clumps may be re- 
garded as physical systems of their own with some of the physics 
necessary to describe them properly not being present in our sim- 
ulations. 



The absolute value of the energy deposition is also resolution 
dependent. It increases by about a factor of 1 .2 if we double the 
resolution. The reason for t h is is l ikely related to the Vishniac 
instability: IVishniac & Rvul dl989l) estimate the wavelength at 
which the growth rate is largest as: 
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where we have plugged in typical values for the column den- 
sity So, the internal pressure P;, and the shell deceleration a. 
This is comparable to our best resolution. Therefore, finer res- 
olution should still trigger strongly unstable Vishniac modes, 
which seem to have an effect on the result. The minimum unsta- 
ble wavelength is predicted to be 0.5/ivi,max- Unfortunately, for 
the present study we did not have the computational resources 
to probe these scales, but this should become possible in future. 
In contrast to the Vishniac instability, the Rayleigh-Taylor insta- 
bility continues to grow faster for smaller wavelengths. Thermal 
conduction would be expected to be important at even smaller 
scales of about 0.01 pc/n, wher e n is the number density in 
the shell dMcKee & Cowidll977l) . Thus, our absolute efficiency 
numbers are lower limits. 

At the level of this accuracy, 3D effects might be important, 
because the shell instabilities should be 3D in nature. Up to the 
first supernova, our simulations are dominated by the wind of 
th e 60 M R star , and m ay thus be compared to the 2D results 
of lFreveretall d2003l) . We find an energy response of at least 
10 per cent, which compares to 9 per cent in the simulation of 
iFrever etaD d2003l) . which is very similar. It might point to some 
effect in the direction that more energy is retained in the ISM in 
3D simulations, but could also be related to numerical details or 
the slightly higher density they use. 

The general shape of the curves is however well converged 
(compare Figure [TT]). As a further check, we have also resimu- 
lated run 3S0 at the resolution of 3Sl-hr. The energy deposition 
ratio between the two high resolution simulations is very sim- 
ilar to the one of the low resolution versions. Additionally, we 
show show below that the relative shapes of the curves are read- 
ily understood. We therefore believe that the relative trends of 
the energy deposition we report here are reliable. 

We find that the Vishniac instability dominates the shell evo- 
lution. We show that the instability in our simulations is con- 
nected to the shells' overdensity and to non-r adial motions in 
the sh ells, i n agreement with the predicti ons of IVishniac & Rvul 
dl989h and iMac Low & Normanl d!993l) . Limiting the shell's 
overdensity by e.g. magnetic fields would therefore directly af- 
fect the Vishniac instability. 

From the column density plots (Figures [3] H and |S), it is 
obvious that the observational appearance of the shell is dom- 
inated by the Vishniac instability: If the shells were smooth, and 
the maximum density would increase with resolution as seen in 
our simulations (Figure [5}, one would expect that the shell gets 
thinner with finer resolution, as the smaller cells allow higher 
compression. Yet, we find a radially averaged shell width of 
tens of per cent of the outer radius independent of resolution 
(FigureQjJi. In the low resolution simulation, much of the width 
is due to the large scale distortion influenced by the grid direc- 
tions, for the high resolution simulation the width is due to small 
wavelength modes. 

In their survey of 322 interstellar bubbles, IChurchwell et al.l 
d2006l) find typical shell widths of 20-40 per cent of the outer 
radius. Thus, it seems unlikely that the development of the 
Vishniac instability is frequently impeded by anything, e.g. lim- 
ited compression due to magnetic fields, as this would again 
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make the shells thin. In other words, in order to study the effects 
of magnetic fields one probably needs much higher numerical 
resolution than adopted in our models. 

The column density should give a rough indication on ob- 
served morphologies. From the corresponding plots, we find 
that the Vishniac instability should also lead to observable fil- 
amentary structure inside the bubbles. This seems to be the case 
for some she lls associated with supernova remnants (e.g. Crab, 
lHesterll2~008l) . which confirms the above analysis. More detailed 
comparison would of course be interesting. 

We find that the best way to inject energy into the ISM, i.e. 
to achieve a high energy response is a continuous, steady en- 
ergy injection. Supernovae dissipate their energy within about 
1 Myr. We show the kinematics for run 3S0 (all stars at same po- 
sition) in Figure [15] After each supernova, the shell accelerates 
significantly. This means more kinetic energy in the shell. Yet 
the increased expansion leads to fast adiabatic pressure loss of 
the shell interior. The increased kinetic energy is quickly dissi- 
pated at the leading radiative bow shock, as long as it is strongly 
supersonic. In contrast, the energy fraction deposited in the ISM 
in the wind phase remains roughly constant at ten per cent. Thus, 
retaining the injected energy in an interstellar bubble requires 
continuous energy injection. 

The energy tracks of merging bubbles are entirely dominated 
by these shell kinematics effects. For example, in run 3S1, the 
merging process has clearly set in at 2 Myrs (compare the high 
resolution version, Figure [3]) and continues for a few Myr there- 
after. Yet, the energy track for this time interval is indistinguish- 
able from run 3S2 (different positions of the stars) and even from 
3S0 (no shell merging because drivers are at same location) and 
the sum of S25, S32 and S60 (no shell merging because the stars 
are sufficiently far away, realised by having them in different 
simulations). 

The lower mass stars inject much less energy into the ISM 
than the high mass star. They produce smaller but still substantial 
bubbles in the ISM. The invariance of the energy tracks within 
the first few Myrs therefore argues for some independence of 
the result of the detailed initial conditions. Some mass loading is 
taking place during shell merging. Yet, also this has no observ- 
able effect on the energy track. 

Exploding a supernova in a superbubble and not in its own 
wind bubble leads to weaker radiative losses: Each supernova 
shock heats first the bubble interior. It then makes a difference 
how large the respective bubble is in communicating the thermal 
energy to the shell: For larger bubbles, the heat energy is dis- 
tributed over a greater volume. Thus the overpressure is smaller. 
The force on the shell is correspondingly smaller. Hence, shell 
acceleration and adiabatic losses of the bubble interior happen 
on a longer timescale. This is the reason for the longer energy 
decay timescale for each subsequent supernova. Consequently, 
after a supernova, the energy decays fastest if the bubbles remain 
isolated, as each star has a small bubble of its own. 

Off-centre explosions are another significant effect for the 
energy tracks: The first supernova always explodes roughly in 
the middle of the superbubble. This must of course be so at 
least for coeval stars, since its parent star also has the highest 
energy output and is the dominant driver of the superbubble be- 
fore it explodes. The energy tracks of the simulations with dif- 
ferent configurations show little difference up to the point when 
the second star explodes. This happens necessarily significantly 
off-centre. The explosion accelerates first and most efficiently 
the parts of the super-shell which are most nearby (compare the 
pressure maps in Figure[T4b. Yet, if the bubbles are fully merged 
at the time of the explosion (3S1) the effect is only at the per 



3S0: shell kinematics 




5 10 15 20 25 30 

Time / Myr 



Fig. 15. Shell kinematics (top: radius, middle: velocity, bottom: 
acceleration), as functions of time for run 3S0. The velocity 
points are averaged over time intervals of varying length, which 
correspond to shell radii differences of at least 2 cells. The shell 
velocity converges towards the ambient sound speed (red line). 
Each supernova leads to a significant acceleration of the shell 
(black crosses, bottom plot), followed by a comparably strong 
deceleration (red stars). 

cent level. This is due to the high sound speed within the bubble, 
which communicates pressure differences quickly. We notice a 
considerable effect on the energy track for run 3S2, where the in- 
dividual bubbles are still well identifiable at the time of the final 
supernova. 

Thus, especially where the shells are not yet fully merged 
at the time of explosion, the off-centre location leads to a cer- 
tain extent to a behaviour closer to the isolated bubble case. 
Therefore, the energy tracks (Figure[T2l of runs 3S1 and 3S2 es- 
sentially do not leave the range spanned by the isolated bubbles 
case (S25+S32+S60) and the cospatial parent star case (3S0). 

Another finding which might seem curious is that all the en- 
ergy tracks in Figure [12] converge at late times. Long after the 
energy injection has ceased, the energy of the affected gas is 
dominated by the kinetic energy of the shell. Because the swept 
up mass is dominated by the action of the 60 M e star and the 
final shell velocity is always similar to the sound speed of the 
ambient medium, the overall energy increase is very similar in 
all simulations. 

Population synthesis of stellar groups/subgroups combined 
with energy inje ction data from stellar evolutionary models 
(IVoss et al.ll2009l) show that the wind energy dominates within 
the first few Myrs after the star formation event. Later, the en- 
ergy input is dominated by supernovae. Obs erved subgroups 
have an age difference of order a few Myrs (IVoss et al.ll2010l 
120121) . Thus, it appears possible that the energy response (com- 
pare equation ([T|i) is kept high for > 10 Myr by the wind con- 
tributions of different subgroups coming in at slightly different 
times. Observat ions find energy responses of about ten per cen t 
or higher (e.g. lOev & Garcfa-Segural l2004t IVoss et all |2012|) . 
This agrees very well with the results in the wind phase of our 
highest resolution run and might suggest that additional effects, 
which are not taken into account in our simulation and which 
we believe should only increase the energy response, may not 
dominate. 

A similar energy response ha s also been inferred observa- 
tionally for galactic winds (e.g. IVeilleux et al.1 120051) . though 
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only the supernova energy has been taken into account for the 
calculation. Galactic winds are thought to arise as a final merg- 
ing stage from central superbubbles in star-forming galaxies. If 
one wants to keep the energy response high in order to match 
the constraints from the galactic wind observations, the individ- 
ual bubbles should merge early in order to have as constant an 
energy input rate as possible. 



5. Conclusions 

We have simulated isolated interstellar bubbles and emerging su- 
perbubbles which form from adjacent interstellar bubbles with 
stellar distances of order tens of pc. Thus, our simulations ap- 
ply, within the limitations outlined in Section |4] above, well 
to hie rarchically cluste red star formation c omplexes like the 
Orion (IVoss et alj|2010l). Sc orpius-Centaurus dDiehl et al.ll2010l) 
or Carina dVoss et al.ll2012l) regions. 

We find in our simulations that up to about the second su- 
pernova the total energy of superbubbles is not strongly depen- 
dent on the spatial configuration of the group of parent stars, 
including zero and infinite distance. Off-centre energy injection 
reduces the ISM energy response significantly only, if the in- 
dividual bubbles are not yet fully merged. Thus, from before 
the second supernova onwards the energy response is higher for 
more closely packed configurations. We find on average about a 
factor of two difference in energy response between the isolated 
stars-case and the cospatial stellar configuration. 

Supernovae increase the ISM energy only for very small 
timescales of about 1 Myr, increasing with the size of the su- 
perbubble at the time of the explosion. After that time, the re- 
tained energy is smaller than immediately before the supernova 
(FigurefTOli. The energy response drops by a factor of two shortly 
after the supernova compared to the main sequence wind phase. 
Our simulations are quite realistic regarding the tim e intervals 
in be tween subsequent supernova events (compare IVoss et all 
l2009h . Thus, we conclude that for realistic star clusters energy 
is build up in the wind phases. Supernovae lead to large short 
term energy variations, but only keep up the bubble energy in 
the long run, at a roughly constant level. 

We also find that supernovae that explode inside larger bub- 
bles have a longer energy decay time. The 60 M Q star has pro- 
duced a bubble of > 80 pc diameter at the time it explodes. Thus 
in order to obtain a physically sound feedba ck model, which is 
curre ntly lacking in studies of disk galaxies ( IScannapieco et al] 
1201 lb . it seems essential to account for the wind phase. Further, 
since the energy deposition does essentially not depend on the 
spatial configuration of the stars, up to stellar distances of about 
30 pc in our simulations, it seems reasonable to use stellar clus- 
ters as fundamental feedback units, not individual stars, or in 
other words superbubbles rather than individual bubbles of indi- 
vidual stars, at least for a clustered star formation mode, which 
should according to our simulations be more efficient for feed- 
back purposes. 

We have verified by comparison to theoretical work that the 
appearance of our wind shells is dominated by the Vishniac in- 
stability, which is here for the first time prominently seen in a 
3D simulation. High resolution is essential to obtain the neces- 
sary shell overdensities which are crucial for the development of 
the instability. This effect widens the shell significantly in col- 
umn density plots, which we suggest may explain the large ob- 
served shell widths of 20 per cent of the outer radii and more. 
It also produces filamentary structure in the shell which is also 
well visible in our column density plots. We conclude that fil- 



amentary structure inside interstellar bubbles may be related to 
the Vishniac instability. 
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